Breathing Current Domains in Globally Coupled Electrochemical Systems: 
A Comparison with a Semiconductor Model 
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Spatio-temporal bifurcations and complex dynamics in globally coupled intrinsically bistable elec- 
trochemical systems with an S-shaped current-voltage characteristic under galvanostatic control are 
studied theoretically on a one-dimensional domain. The results are compared with the dynamics 
and the bifurcation scenarios occurring in a closely related model which describes pattern forma- 
tion in semiconductors. Under galvanostatic control both systems are unstable with respect to the 
formation of stationary large amplitude current domains. The current domains as well as the homo- 
geneous steady state exhibit oscillatory instabilities for slow dynamics of the potential drop across 
the double layer, or across the semiconductor device, respectively. The interplay of the different 
instabilities leads to complex spatio-temporal behavior. We find breathing current domains and 
chaotic spatio-temporal dynamics in the electrochemical system. Comparing these findings with the 
results obtained earlier for the semiconductor system, we outline bifurcation scenarios leading to 
complex dynamics in globally coupled bistable systems with subcritical spatial bifurcations. 

PACS numbers: 82.40.-g, 82.45.-h, 72.20.Ht, 05.70. L 



I. INTRODUCTION 

The focus of research in nonlinear dynamics has 
evolved from temporal instabilities over simple spatial 
patterns to complex spatio-temporal behavior and the 
control or synchronization of such dynamics. Complex 
spatio-temporal behavior in reaction-diffusion equations, 
which is in a wider sense the class of equations dealt 
with also in electrochemistry, might be found when in- 
stabilities breaking time and space symmetries inter- 
act. A generic case is the interaction of Turing M and 
Hopf bifurcation in a two component activator-inhibitor 
system in which the involved species diffuse. Com- 
plex spatio-temporal dynamics has been found near this 
codimension-two point theoretically ||, |] as well as 
experimentally & @, [7]] . 

In electrochemical systems that can be described by a 
two component model one variable typically is of elec- 
trical nature and the associated transport mechanism is 
migration rather than diffusion ^| . The decisive vari- 
able for the dynamics of the electric circuit is the double 
layer potential 4>dl, measuring the voltage drop across 
the interface between the working electrode and the elec- 
trolyte solution ]10| . Local perturbations in the double 
layer potential are mediated through the electric field in 
the electrolyte. Thus, spatial inhomogeneities in the dou- 
ble layer potential are felt not only by its nearest neigh- 
bors, but by a whole range of neighboring sites which 
makes the coupling nonlocal pf. The degree of non- 
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locality depends on the geometry of the electrochemical 
cell, most importantly on the positions of the working 
(WE), counter (CE) and reference electrode (RE) with 
respect to each other [Q. Furthermore it can be shown 
that the galvanostatic operation mode (constant current 
control) introduces an additional global coupling into the 
system [fl3|| . 

The role of the second variable in two component 
activator-inhibitor systems in electrochemistry is played, 
e. g., by the chemical concentration of the reacting 
species in the double layer or by the density of adsorbed 
molecules on the WE. 

Over the last decade global coupling has been an active 
area of research. Global coupling is present in systems 
that are subject to external control, e.g. via an electric 
circuit (such as in electrochemical 0, OL OL j^, |U| 17 , 

H H, & semiconductor © H M M E§ M- lj- £$ 

and gas discharge |2!| systems) or via the electric control 
of the temperature in catalytic reactors |50|, [nj ||[ |||, 
p4| , p5[ . But global coupling may also be due to transport 
processes that happen on time scales much faster than all 
other relevant time scales in the system, e. g. fast mixing 
in the gas phase ||, |7|, || [|| goj, fllj. A variety of 
other systems are described by dynamics that include 
global coupling, e.g., ferromagnetic [Q, biological J43| , 
and chemical systems in which the global coupling can 
be light induced 44, 4a|. Abstract theoretical models are 
discussed, e. g., in |§. 

Results regarding electrochemical systems with global 
coupling have been reported for systems with an N- 
shaped current-voltage characteristic (termed N-NDR 
systems: N-shaped negative differential resistance) for 
different types of global coupling. In these systems the 
double layer potential acts as an activator and global cou- 
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pling introduced by the galvanostatic control mode was 
shown to accelerate front motion Jll], [l4| thus having 
a synchronizing effect on the spatial dynamics. Desyn- 
chronizing global coupling of the activator was shown to 
stabilize potential fronts, leading to two stationary po- 
tential domains Jl5|, ^9). Also the formation of pulses 
and standing waves was observed |l^, |5(| . 

In electrochemical systems with an S-shaped current- 
voltage characteristic (S-NDR) the roles of activator and 
inhibitor are reversed, leading to a global coupling of the 
inhibitor under galvanostatic conditions. This leads to 
the opposite effects opposed to N-NDR systems, i.e. cur- 
rent domains that are stabilized by the constant current 
constraint 0. Similar results on accelerated and de- 
celerated fronts in globally coupled semiconductors with 
S- or Z-shaped current-voltage characteristics have also 
been obtained Q. 

In the present paper we focus on the latter case of 
global coupling of the inhibitor in an electrochemical sys- 
tem with an S-shaped current-voltage curve. Further- 
more, we consider systems with high electrolyte conduc- 
tivity. In such systems the migration coupling is so effi- 
cient that any spatial variation in </>dl can be neglected, 
which results in an additional global coupling Eh]. The 
set of equations to be investigated is thus of the general 
form: 



T <fc 



at 



<?(</>DL, (9)g) 



(1) 
(2) 



where 9 stands for the activator variable, whose dynamics 
comprises an autocatalytic chemical step. The angular 
brackets denote the spatial average over the spatial do- 
main G. f is autocatalytic in 9; g exhibits a monotonic 
characteristic as a function of 0dl and 9. T^ DL (e) de- 
note the characteristic times for changes in </>dl and 9, 
respectively. 

A formally very similar set of equations describes 
the dynamics in bistable semiconductor systems oper- 
ated via an external load resistance j2l], |22|, ^8], j5l| . 
The formation and dynamics of current density pat- 
terns in bistable semiconductors was extensively stud- 
ied (2^, p(| [32| |53|, In this respective class of semi- 
conductor systems the current-voltage characteristic also 
resembles the shape of an S, which points to the fact 
that the roles of the dynamic variables are very similar 
to the electrochemical model: The voltage drop u across 
a semiconductor device acts effectively as inhibitor, and 
it is subject to a global constraint imposed by the ex- 
ternal electric circuit. The role of the activator variable 
might be played by different physical quantities, such as 
the electron temperature [El[, the concentration of ex- 
cess carriers [E5L t he charge density in resonant tunnel- 
ing structures |54 [56| , |5?j (note that for bistable reso- 
nant tunneling structures the current-voltage character- 
istic is Z-shaped resulting in an activatory, not inhibitory 
effect of the global constraint), the voltage drop across 



pn-junctions in thyristors |25|, Q or the interface charge 
density in a heterostructure hot electron diode (HHED) 
p3[ . The dynamic equations are of the same structural 
form as eq. (|l]),(^); only the local nonlinear functions / 
and g differ from the electrochemical model. 

For the current density dynamics in a class of mod- 
els originally derived for the HHED in one or two spa- 
tial dimensions under galvanostatic (current-controlled) 
conditions, interesting complex spatio-temporal patterns 
termed "spiking" and "breathing" current filaments were 
found . Recently, a sufficient condition for the on- 

set of such complex spatio-temporal dynamics was given 
& 

Realizing the obvious similarities, we show in this pa- 
per that the methods (e. g. for analyzing the dynamics) 
developed for the semiconductor system can be applied 
to gain insight into the interaction of different instabili- 
ties in the electrochemical system. Results regarding the 
possibility of the occurrence of complex spatio-temporal 
behavior and the mechanisms that lead to such behav- 
ior are given. It is emphasized whether the different 
dynamical regimes depend upon the general structural 
form of the equations, especially regarding the influence 
of global coupling, or if they are due to special proper- 
ties of the underlying physical or chemical system, and 
thus the local dynamics. Hence a comparison of electro- 
chemical and semiconductor systems gives considerable 
insight into generic complex dynamics of globally cou- 
pled bistable systems. 

The paper is organized as follows. In the next section 
we introduce the electrochemical model, discuss its im- 
portant parameters and the mechanisms leading to global 
coupling in the model. In the following section we char- 
acterize the dynamics of the model by linear stability 
analysis along the lines developed for the semiconductor 
model and by numerical simulations. In the discussion we 
compare the important features of the two models. The 
mechanism leading to complex spatio-temporal behavior 
in both models is different and this difference is explored 
in this section in some depth. We summarize our results 
in the last section and give a short outlook to applications 
in terms of experimental verifications and transfer of the 
electrochemical results to the semiconductor model. 



II. MODEL 

The central variable in electrochemical pattern forma- 
tion is the double layer potential (/>dl , the potential drop 
across the interface between the WE and electrolyte so- 
lution. The dynamic evolution equation for 0dl can 
be deduced from the local charge balance at the elec- 
trode/electrolyte interface. 

To make things as transparent as possible and to fa- 
cilitate later comparison with the semiconductor model, 
we employ the idealized geometry shown in fig. § WE 
and CE are equally sized rectangular plates positioned 
parallel to each other in a box-like cell with otherwise in- 
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FIG. 1: Schematic setup of the electrochemical system. A 
constant current io is applied in a electrochemical cell con- 
sisting of a rectangular working electrode (WE), electrolyte, 
and a rectangular counter electrode (CE). WE and CE form 
top and bottom of the box-like cell with otherwise insulating 
walls. (j>uL is the voltage drop across the interface. i t and 
ic symbolize reaction current density and capacitive current 
density, respectively. 



undergoes a first order phase transition due to lateral in- 
teractions of the adsorbate molecules 59, S0|. The trans- 
formations leading to dimensionless units differ from the 
ones given in |59f| ; the derivation is given in Appendix |a|. 
Note the non-polynomial nature of the function /. 
The dimensionless set of equations is thus 



DDL 



dt 
86 



= 7 (io-(l- (0))e^) 



5 -"("-* 



-«j(0,0dl) 



p9e 



iu(0,0dl) 



(6) 

d 2 e 

dx 2 
(7) 



with 



w(0, 0dl) = v<j>\ 



DL 



sulating walls. This geometry imposes no-flux boundary 
conditions for 0dl and 6; there will be no spatial inho- 
mogeneities of the electric field at the interface imposed 
by this geometry. 

For very high electrolyte conductivities a, spatial in- 
homogeneities in the double layer potential are damped 
out very fast via the efficient coupling through migration 
currents. It follows that spatial variations of 0dl can be 
neglected. This effectively introduces a global coupling 
in the system, since local perturbations in c/>dl are felt 
instantaneously in the whole double layer. 

In the following we additionally assume current con- 
trolled conditions. Galvanostatic control is known to 
introduce an additional global coupling into the system 
Assuming a specific double layer capacitance C, 
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the dynamic equation for the double layer potential reads 



C- 



dt 



= ~«r(0DL, (0)) + io, 



(3) 



where i r (0DL, &) is the reaction current density and io de- 
notes the imposed current density. The activator variable 
8 describes the evolution of the coverage of the WE by an 
adsorbate or the concentration of a chemical species in 
the reaction plane. Its dynamics will be modeled by an 
equation of the form (g) , where we restrict our system to 
one spatial dimension (Id) since the qualitative behavior 
should also be captured on Id domains. Id domains also 
resemble the situation of a very large aspect ratio of the 
rectangular domain, where one spatial dimension is too 
small to allow for spatial instabilities and can thus be 
eliminated. 

We use the following model functions for the local dy- 
namics of the activator and the reaction current density 



!>DL, 



= (1 



(4) 
(5) 



originally derived to describe pattern formation observed 
in a reaction, in which a reaction inhibiting adsorbate 



subject to the boundary conditions 



dx 



= 0. 



X—Q.7T 



We normalize space to the interval [0, tt] for computa- 
tional convenience and thus 



(6>) = i [ 9(x) dx. 
7T Jo 



This leads to the proportionality of the parameters 
/i,7 oc L 2 ; fi and 7 still can be changed independently 
since also other physical quantities enter these parame- 
ters (cf. Appendix |A|). 

In fig. (|2|a) the nullclines of the system are depicted for 
a current density io that is set in the range of the neg- 
ative differential resistance in the current-voltage char- 
acteristic (see fig. ||b). The S-shaped current- voltage- 
characteristic is depicted together with the load line 
i = io in fig. |^b. This physically more intuitive (i-(puh)- 
plane representation will be used in the following. 

The parameters v, p and g are fixed throughout this 
paper at the values v = 0.025, p = 0.5 and g = —2.4 
(cf. Appendix ^). The dynamics is determined by the 
model parameters //, essentially proportional to L 2 , the 
relaxation time ratio of activator and inhibitor 7/ /i (inde- 
pendent of L) , and the general excitation level controlled 
by the imposed current density io- The relaxation time 
ratio can be accessed easily via the concentrations of the 
reacting and adsorbing species; io is set by the galvano- 
static control unit. 



The numerical results discussed in section III B were 



obtained using pseudo spectral decomposition in space 
[ ]6l"[ employing 15 spatial cosine- modes (the results do 
not change when a larger number of modes is chosen). 
For the integration in time the routine lsode p2] | and for 
continuation of steady states and limit cycles the package 
AUTO ||S3| was used. 
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TABLE I: Abbreviations for bifurcation points 

h Hopf bifurcation of the homogeneous steady state 

d domain bifurcation of the homogeneous steady 

state 

sn-d saddle-node bifurcation of domains 

hd Hopf bifurcation of the domain state 

snp saddle-node bifurcation of breathing domains, i. 

e. periodic orbits 
DH domain-Hopf codimension-two point (d and h) 

TB Takens-Bogdanov codimension-two point (sn-d 

and hd) 

DHD degenerate Hopf bifurcation of domains (snp and 
hd) 



(c) 




FIG. 2: (a) Nullclines of the model (§),(§) for an imposed 
current density in the autocatalytic regime (solid line 9 = 0, 
dashed line <^dl = 0, io = 1300, for the other parameters see 
Appendix |a|). (b),(c) S-shaped current- voltage curve together 
with the load line io (jo) for the electrochemical (eq. @,(Q)) 
and the semiconductor (eq. (^3|),(|l^)) system, respectively. 



Subscripts denote partial derivatives with respect to 
the subscripted variable and evaluation at the steady 



state (e. g. f e 



ill 

,w 



). For brevity we denote 



(«>E 



The stability of the fixed point with respect to homo- 
geneous fluctuations (n = 0) can be determined by in- 
specting 



det J 



fe 



d6 s 



DL 



\ a<?DL 

dl r (0 ss (^DL),0DL) 



-7M/< 



III. STABILITY ANALYSIS AND 
SIMULATIONS 

A. Homogeneous Steady State 

In this section we consider the spatially uniform 
fixed points of the system @,(0) and their bifurca- 
tions. The uniform steady state (0q L ,6* ss ) is given 
by i r ((/)^ L ,9 ss ) = i , /(0dl;^ ss ) = an d corresponds to 
the homogeneous S-shaped current- voltage characteristic 
(fig. ||b). Perturbing the steady state with a perturbation 
(<5</>DLe At , 59 cos(nx)e xt ) (consistent with the boundary 
conditions), the temporal evolution of the perturbation 
is given by the eigenvalues of the Jacobian matrix J 

trJ , /(trJ) 2 " " 
Ai,2 = — ±y^--detJ 

and stability (ReA < 0) implies that (det J > A 
trJ < 0). The Jacobian reads 

j = ( ~Wt ~lh e \ 



and 

trJ = fifg - ~fcr r . 

Obviously det J > in general since [i, 7 > and 
/ /' r(rS it DL) ' 0DL) < 0, which follows from the fact 
that the branch of negative differential resistance 

( d ' r(9 l0DL ),0DL) < 0) is caused solel y b y the activator 
variable 6, equivalent to saying that er r > in general. 

However, tr J might change sign on the NDR-branch 
since fe>0 and er r > 0, which leads to an oscillatory 
instability (Hopf bifurcation) of the homogeneous steady 
state (denoted by a superscript "h" , cf. table |) at 




Thus for 7/// < (7/yLt) (low concentration of the reacting 
species or high concentration of the adsorbate) the homo- 
geneous steady state is unstable in a certain io-interval, 
since fe<y~ x depends on the imposed current density via 
the steady state condition. When plotting the critical 
value ("i / iif 1 as a function of the imposed current den- 
sity fig. |^a is obtained. 
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parameter value is thus: 



0.8 1.0 1.2 1.4 1.6 

10" 3 x i 

(b) 50 p 



25 



H m »=3.54 



0.8 1.0 1.2 1.4 

10" 3 x i„ 



1.6 




FIG. 3: (a) Location of the Hopf-bifurcation of the homo- 
geneous steady state in the (7//i-io)-parameter-plane for the 
electrochemical model (^|),(0)- F° r t/A 4 > 2.2 • 10 -4 the sys- 
tem is stable with respect to homogeneous fluctuations, (b) 
Threshold for the spatial instability of the uniform steady 
state in the (/x-io)-plane. For system sizes smaller than 
Mmin = 3.54 the system is stable with respect to spatial fluc- 
tuations, (c) Critical system size L„. of the spatial instability 
for the semiconductor model (eq. (|13|),(p^|)) as a function of 
the imposed current density jo- 



For 7//1 > (7/^)Jnax = 2-2 ■ 10~ 4 there are no oscilla- 
tory solutions for any %q and for 7//1-C (7/A*)m a x the 
oscillatory instability takes place close to the turning 
points of the current-voltage characteristic at io = 889 
and io — 1587. 

To determine the stability with respect to spatially in- 
homogeneous fluctuations, it is sufficient to consider the 
activator variable 9, since sinusoidal perturbations do not 
affect the average value of and thus the double layer 
dynamics. Therefore the steady state becomes unstable 
with respect to the n th -mode for 



A* > 



fo 



and the first mode to become unstable is always the mode 
with wavenumber one pif . The wavelength of the first 
unstable mode depends on the system size and is equal to 
2L for Neumann boundary conditions. In the following 
we term this instability domain bifurcation. The critical 



/ = fo 1 - 



(9) 



This critical value is depicted in fig. ||b as a function of 
io. For systems sizes \i < fi m - m = 3.54 the spatial insta- 
bility is suppressed; this defines a natural length scale 
for the system. For system sizes much larger than this 
natural length scale the spatial instabilities occur once 
again close to the turning points of the current-voltage 
characteristic. 

The spatial and oscillatory instabilities may coincide 
in a codimension-two point (Domain-Hopf bifurcation, 
"DH" , cf. table @) if 

7 1 . (10) 



,,DH 



7 



The respective imposed current density value «o (m) is 
defined as the solution of (0) with respect to io- 



B. Homogeneous Limit Cycle and Stationary 
Domains 



In this section we complete the picture of the different 
basic attractors of the model by including limit cycles 
and stationary current domains into our stability anal- 
ysis. Analytical methods fail in most cases since the 
involved bifurcations are either subcritical and thus do 
not allow for an amplitude equation analysis and/or the 
considered system sizes are intermediate, which excludes 
methods like singular perturbation theory [ |65| to describe 
domain interface dynamics. 

For common concentrations and system sizes the dou- 
ble layer dynamics will be much faster than the dynamics 
of the activator. For these conditions the parameters 7 
and \x will be of the order 10 and 100, respectively. It 
follows that in most cases oscillatory instabilities are not 
present in the system and the only nontrivial mode is a 
stationary current domain as depicted in fig. |^a for two 
values of io- This current domain is the final state of the 
system in the spatially unstable regime and the mecha- 
nism leading to such a stationary domain is well known 
(e.g. 0,0): 

The activator is bistable as a function of the double 
layer potential. An overcritical local fluctuation in a 
system without global coupling that is prepared in the 
metastable state would lead to the formation of a tran- 
sition front to the globally stable state. The global con- 
straint, however, forces the system to maintain an aver- 
age current. The system meets this constraint by taking 
on an inhomogcncous state in which two phases coex- 
ist. With other words, the front velocity becomes zero. 
The final state of the system is described by a Maxwell 
type construction: the intermediate, equistability double 
layer potential </>p q L , which is established in the stationary 
structure is determined by the equal-areas rule |2q, p4] 



/(^ eq 



DL> ' 



I d6 = 0. 



(11) 
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FIG. 4: (a) Stable domains for two values of the imposed cur- 
rent density io for the electrochemical model (fi = 25, 7 = 10). 
(b) Bifurcation diagram for fi = 25 and 7 = 10. Shown is </>dl 
as a function of the bifurcation parameter io in the familiar 
current-potential plane. The branch of negative differential 
resistance is unstable (thin dashed line) with respect to do- 
main formation. The domain branches (thick lines) bifurcate 
subcritically (d) near the turning points of the current-voltage 
characteristic. The stable and unstable domain branches 
(solid resp. dashed thick lines) are born in a saddle-node bifur- 
cation of domains (sn-d) . The domain branch can be approx- 
imated by an equal-areas rule, eq. (|ll|), in a huge io interval 
yielding an equistability potential <^ L . 



In fig. the bifurcation diagram with respect to io is 
shown for /1 = 25, 7 = 10. Even though the system size 
is comparable to the interface width, as can be seen 
111 fig. |a the above construction holds for a wide i§- 
interval. However, since the arguments given above ap- 
ply strictly only for infinite systems, deviations near the 
turning points of the current- voltage characteristic of the 
domain are clearly visible. These deviations represent a 
boundary effect. 

States with several domains are unstable due to the 
winner takes all principle [ pl| p?| . Domains with an ex- 
tremum not located at the boundaries are unstable with 
respect to translation and will be attracted by the bound- 
ary. 

Fig. |^b also shows that spatially patterned solutions 
typically bifurcate subcritically from the homogeneous 
state and meet the stable domain-branch in a saddle-node 
bifurcation (sn-d, cf. table Q). The domains remain sta- 
ble in the whole zo-interval in which the current-voltage 
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4 , io = 1000). Shown is the reaction current den- 
^ DL as a function of space and time. At these 



FIG. 5: (a) Bifurcation diagram for the electrochemical model 
(jit = 25, 7 = 3- 10~ 3 ). Apart from the branches shown in 
fig. ^b, a branch of unstable homogeneous oscillatory so- 
lutions (open circles) bifurcates supercritically (h) near the 
turning points of the current-voltage characteristic. Shown 
is the maximum value of </>dl during one oscillation cycle. 
After stabilization through a pitchfork bifurcation, the sta- 
ble homogeneous oscillations (full circles) are separated from 
the stable domains by an unstable inhomogeneous limit cycle 
(open triangles), (b) Typical scenario of an oscillatory insta- 
bility of a domain for lower values of 7 than in (a) (/1 = 25 
7 = 1 • 10 
sity i r = (1 — V)e 
parameter values the oscillatory instability of the domain is 
subcritical and the system finally settles down to homoge- 
neous relaxation oscillations (standard scenario). 



characteristic of the domain exhibits a negative differen- 
tial resistance for these parameter values. This behavior 
can also be rationalized analytically |2^, . The domain 
bifurcation is supercritical only in a small /i-interval close 
to the minimal system size /i m in- 

When fi is fixed at a value [i > /i m in and the double 
layer dynamics is slowed down to 7 below M(7/Ai)m ax , 
the additional mode of homogeneous oscillations becomes 
present in the system. For 7 < ^(7 / /i)^ lax it bifurcates 
supercritically from the spatially unstable state, therefore 
small amplitude oscillations will be unstable with respect 
to spatial fluctuations for any io. With increasing oscil- 
lation amplitude (decreasing 7) the oscillations become 
stabilized in a pitchfork bifurcation (cf. fig. ||a). This re- 
sults in bistability of stationary domains and an uniform 
limit cycle in an intermediate i -'mteTval. The basins of 
attraction are separated by an unstable inhomogeneous 
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limit cycle. 

If 7 is lowered even further, the stationary current do- 
main will become unstable also. This can be rationalized 
by recalling that the stabilization mechanism of the do- 
mains is the positive global coupling on 0dl ■ If the delay 
of the double layer dynamics becomes too large, </>dl can 
no longer control the interface stability. We denote the 
critical value of this oscillatory instability of the domain 
by 7 hd (/i, «o) (cf. table [jj). Numerical simulations show 
that the threshold for an oscillatory instability of the 
current domain lies typically below the threshold for the 
Hopf bifurcation of the homogeneous steady state 

7 hd (AMo) <7 h (AMo). (12) 

This can be understood in the frame of the eigenmodcs 
of the current domain for large system sizes if we recall 
that in absence of global coupling the domain state has 
only one positive eigenvalue that tends to zero with in- 
creasing system size. The respective arguments are given 
i n p3| . The numerical investigations show that relation 
(|l2|)in general holds for small and intermediate system 
sizes also. 

It follows that, in general, the homogeneous relaxation 
oscillations represent an attractor when the domain loses 
stability. The oscillatory instability of the domain is usu- 
ally subcritical; a state close to the domain is eventually 
attracted by the stable homogeneous limit cycle (see fig. 

. This can be regarded as the standard scenario (i. e. 
it exists in a wide parameter range) of a domain instabil- 
ity in globally coupled electrochemical systems with an 
S-shaped current-voltage characteristic. In this case no 
complex spatio-temporal behavior arises in the model. 

C. Breathing Current Domains 

We would expect complex spatio-temporal behavior if 
the branch of inhomogeneous limit cycle solutions that 
bifurcates from the domain- branch at the point of the os- 
cillatory instability of the domain becomes stabilized or 
bifurcates supercritically. In this case the system would 
exhibit bistability between a stable inhomogeneous limit 
cycle and a stable homogeneous one. We did indeed find 
such a situation in the model for comparatively small 
system sizes (fi ~ 10) and relaxation times well below 
the onset of homogeneous oscillations (7//i^7-10 -5 ). 
The instability leading to such complex spatio-temporal 
behavior is shown in fig. ^|. In fig. [f]a the correspond- 
ing bifurcation diagram for /i = 10 and 7 = 7- 10~ 4 is 
depicted. 

Decreasing the imposed current density from values in 
the regime of bistability between a stable domain and 
homogeneous oscillations, the domain branch exhibits an 
oscillatory instability (hd) . The branch of solutions that 
bifurcates subcritically is stabilized via a saddle-node bi- 
furcation of oscillatory domains, i. e. periodic orbits (snp, 
cf . table ||) , which can be seen in the enlarged bifurcation 
diagram, fig. p]b. The spatio-temporal behavior becomes 
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FIG. 6: Oscillatory instability of a domain leading to sta- 
ble, periodically breathing current domains for the electro- 
chemical model = 10, 7 = 7- 1CT 4 , i = 1000). (In this 
simulation a stable domain was prepared, 7 was lowered to 
7 = 7- 10 -4 , and a small random fluctuation was added.) 



TABLE II: Dynamic regimes indicated in fig. || 

(1) One stable homogeneous fixed point. 

(2) Bistability between stable domain and homogeneous 
fixed point. 

(3) Only one stable domain. 

(4) One stable homogeneous limit cycle. 

(5) Stable or unstable homogeneous limit cycle (cf. Ap- 
pendix B) and stable domain. 

(6) Stable homogeneous limit cycle and stable domain. 

(7) Stable breathing current domains (periodic or chaotic) 
and a stable homogeneous limit cycle. 

(8) The Hopf bifurcation of the domain is subcritical, thus 
only stable homogeneous oscillations are present. 

(9) Region in which three attractors exist (cf. fig. for io = 
1010): Stable domains, stable breathing domains, and 
stable homogeneous limit cycle. 



more involved as the imposed current density io is de- 
creased. The limit cycle undergoes a period doubling 
cascade leading to stable chaotic spatio-temporal motion 
(fig. ||). Decreasing io further, a reversed period dou- 
bling cascade occurs which leads again to stable period 
one breathing domains. This branch then ends in a su- 
percritical Hopf bifurcation of the domain very close to 
the saddle-node bifurcation, in which stable and unsta- 
ble domains originate (sn-d). It is interesting to note 
that the dynamic nature of the invariant set that sepa- 
rates the basins of attraction of the two limit cycles is 
changing with increasing imposed current density from 
the unstable stationary domain (saddle point) to an un- 
stable inhomogeneous limit cycle (see fig. 0b for low io). 

The region in the (io-7)-parametcr-planc in which such 
complex spatio-temporal dynamics is found is depicted 
in fig. H for /i = 10. The lines of the Hopf bifurcation 
and the domain bifurcation of the homogeneous steady 
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FIG. 7: (a) Basis bifurcation diagram for stable pe- 
riodic breathing for the electrochemical model (fi = 10, 
7 = 7- 10~ 4 ). The oscillatory branch of the homogeneous 
limit cycles bifurcates supercritically before the spatial in- 
stability and thus homogeneous oscillations are stable nearly 
in the whole io-interval (full circles). The equal areas rule, 
eq. (|ll|), fails for this system size. The domain branch (thick 
line) is unstable in a region of negative differential resistance 
(dashed thick line) near the lower saddle node bifurcation of 
domains. Marked with open triangles is an unstable inhomo- 
geneous limit cycle. It is born in a pitchfork bifurcation of 
the homogeneous limit cycle at high io and terminates in the 
unstable domain branch, (b) Enlargement of the bifurcation 
diagram at the lower turning point. Here also the branches 
of the inhomogeneous breathing mode are shown (diamonds) . 
The breathing mode bifurcates subcritically (hd) from the do- 
main branch at higher io (open diamonds) and stable breath- 
ing (full diamonds) originates in an snp. In the projection 
of the limit cycle on the double layer potential it gets close 
to the homogeneous steady state but not in real phase space 
(cf. text). In the current density interval between approx. 
io = 1000 and io = 975 the inhomogeneous limit cycle under- 
goes a period doubling cascade leading to chaotic breathing 
(open diamonds). 
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FIG. 8: (a) Periodically breathing domains with period two 
at io = 990 and (b) chaotically breathing domains at io = 980 
Qu = 10, 7 = 7-10- 4 ). 
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FIG. 9: Existence region of stable periodic breathing in the 
(io-7)-control-parameter plane (hatched region) for the elec- 
trochemical model (fj, = 10). The main dynamic regimes, 
characterized by attractors, are indicated by the numbers 1— 
9; the attractors are given in table Shown in solid lines 
are the lines of the domain bifurcation (d) and the saddle- 
node bifurcation in which the domains originate (sn-d) (both 
independent of 7). The domain bifurcation and the Hopf 
bifurcation of the homogenous steady state (h, solid) inter- 
sect in a Turing-Hopf-type codimension-two point (DH). The 
dashed line shows the oscillatory instability of the domain 
(hd). Denoted by snp (solid line) is the line of the saddle 
node bifurcation of periodic orbits, i. e. breathing domains. 



state and their intersection point (DH) are shown. The 
main regions that were discussed above (and in part also 
exist for different values of fi) are indicated. Note the 
existence of three codimension-two points: The point in 
which domain an d Hop f bifurcation coincide (DH) was 

At the DH the system has a 



discussed in sec. [II A 



pair of purely imaginary eigenvalues and a real eigen- 



value equal to zero |6(|. Unfoldings of the DH have a 
further fine structure as discussed in Appendix it is 
not shown in fig. ^| for clarity. Denoted by TB is the point 
where saddle-node and Hopf-bifurcation meet (Takens- 
Bogdanov point) |56|. Note that in our case both bi- 
furcations involve inhomogeneous steady states (i. e. do- 
mains) rather than homogeneous solutions. Left of the 
TB two saddle fixed points with one and three unstable 
directions, respectively, originate from the saddle-node 
bifurcation, right of it a saddle fixed point and a stable 
node. Again the fine structure, most remarkably a homo- 
clinic bifurcation that should be present in the vicinity 
of the TB, is omitted. The third codimcnsion-two point 
is a degenerate Hopf bifurcation of domains (DHD), in 
which the saddle-node bifurcation of periodic orbits (snp) 
coincides with the Hopf-bifurcation of the domain (hd). 

We omitted in the bifurcation diagram (fig. |fj) some 
of the branches mentioned above. Furthermore there are 
indications of the presence of additional bifurcations that 
determine the exact location of the lower boundary of the 
regime of complex behavior. Its detailed study is beyond 
the scope of this paper. 



IV. COMPARISON AND DISCUSSION 

In this section we compare the different dynamic in- 
stabilities and regimes described in the previous sec- 
tion with results obtained earlier for the semiconductor 
model. The semiconductor model used has the (nondi- 
mensionalized) form 

u — a d 2 a 
a=- - 0.05a +— 13 

(u - a) 2 + 1 dx 2 y ' 

ii = a(j Q -u+ (o)l), (14) 

where u denotes the potential drop across the semicon- 
ductor device (corresponding to 0dl) and a describes the 
interface charge density in the HHED (corresponding to 
6). The system length is L and thus (a)^ = L^ 1 J Q adx. 
The current-voltage characteristic of the HHED is given 
by j = u — a. It also has the shape of an 'S' (fig. fy). If 
space is rescaled to the interval [0, it], the model exhibits 
the same structural dependence on three parameters as 
eqs- ©,(0) 

\ / u — a \ d 2 a 

\(u-a) 2 + l J dx 2 y ' 

« = 7 (s) (io-«+(aM, (16) 

with /j,( s ) = (^-) 2 and 7W = (^-) 2 a. These parameters 
can be interpreted in the same way as in the electro- 
chemical model. 

The two models possess equivalent basic modes: The 
branch of negative differential conductivity is unstable 
with respect to spatial perturbations for sufficiently large 
system sizes L > L m j n (cf. fig. |^c) and to homogeneous 
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FIG. 10: (a) Spiking current filament in Id (L = 40, 
a = 0.035, jo — 1.2). (b) Bifurcation diagram for the semi- 
conductor system for complex spatio-temporal dynamics 
(L — 40, a = 0.06). Shown is the potential drop across the 
semiconductor u resp. the maximum u during one oscillation 
versus the imposed current density jo at the lower turning 
point of the current-voltage characteristic. In the current- 
interval shown by the dashed lines no trivial state of the sys- 
tem is stable. The lower boundary is the spatial instability 
of the homogeneous steady state (thin lines) and the upper 
one is the oscillatory instability of the filament (thick lines). 
Homogeneous oscillations are not present in this current den- 
sity interval, they bifurcate at higher current density values 
(open circles in the upper right corner). The resulting in- 
homogeneous oscillations (diamonds) that bifurcate subcriti- 
cally from the stable domain branch are born by a saddle-node 
bifurcation of periodic orbits. 



oscillations for sufficiently slow dynamics of the voltage 
drop u (small a). However the temporal instability of 
the filament may lead to qualitatively different spatio- 
temporal dynamics. Apart from the breath ing mode that 
the semiconductor models also exhibit, [^6[ |53|, |67j the 
system displays a complex spatio-temporal mode termed 
spiking (see fig. |lO|a) [ pjl |52| . 

This mode evolves because the spatially inhomoge- 
neous limit cycle that constitutes breathing comes even- 
tually, with decreasing 7, very close to the homogeneous 
fixed point. This points to a structurally different dy- 
namic regime as compared to the electrochemical system 
and facilitates the formulation of a sufficient condition 
for the occurrence of complex spatio-temporal dynam- 
ics ^(|. In the following this will be explained in some 
detail. 

Consider the bifurcation diagram of the semiconduc- 
tor model for parameter values at which complex spatio- 
temporal dynamics is found (fig. [To|b). Let us de- 
note by join), Jo(m>7) an d Jo d (M,7) the parameter val- 
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ues at which the spatial instability of the homogeneous 
steady state, the oscillatory instability of the homoge- 
neous steady state and the oscillatory instability of the 
filament, respectively, occur. For an interval of imposed 
current densities jo n o trivial state is stable, since, in con- 
trast to the electrochemical model, homogeneous oscilla- 
tions are not present in the system for imposed current 
density values within this interval, however the filament 
is already oscillatorily unstable (j hd (^j,,io) > 7 h (/^,io))- 
Thus a sufficient condition for complex dynamics is 

j h (/i,7)>j d M A Jo hd (M,7)>Jo d (M)- 



The limit case j hd (/., 7) = ftfa 7) = J$(v) = J DH (^) 
can be reformulated as a condition for the timescale of 
the inhibitor 7 such that the condition can be tested for 
different system sizes p6|: 



7 hd ^))>7 Ui V)- 



(17) 



The above inequality becomes clear, if one considers that 
the oscillatory instability of the filament is shifted toward 
higher imposed current density values when lowering 7, 
whereas the Hopf bifurcation point of the homogeneous 
steady state behaves in the opposite way and the spatial 
instability does not depend on 7. 

In fig. 11 both critical timescales are plotted for 
both models. For the electrochemical model the criti- 
cal timescales are also shown for the upper part of the 
S-shaped current-voltage characteristic. The above ar- 
guments do equally apply for this region. As indicated 
by the hatched region for the semiconductor system in 
fig. [Tl|a, condition ( |l7| ) is fulfilled for a large interval 
of system sizes L (resp. /i^) for the lower part of 
the S-shaped current-voltage characteristic. Apart from 
spiking, a broad variety of periodic and chaotic spatio- 
temporal modes has been found in this interval 1 26 . Con- 
dition jTfl) is never found to hold for the upper part for 
the semiconductor model (not shown). It can be seen in 
fig. |ll] b and [TT| c for the upper and lower part of the S- 
shaped current-voltage characteristic, respectively, that 
condition (|l7|) is apparently never fulfilled in the electro- 
chemical system for any system size. 

Thus also the absence of spiking in the electrochemi- 
cal system is easily understood; spiking evolves when the 
breathing mode eventually comes very close to the plane 
of homogeneous modes which constitutes a stable focus 
in this plane. The relaxation close to the homogeneous 
fixed point in the plane of the homogeneous modes leads 
to the small, almost homogeneous, oscillations and then 
the spike evolves again as the trajectory leaves the plane 
of homogeneous dynamics along the unstable direction 
of the homogeneous fixed point (cf. fig. 



10). In the elec- 



trochemical system the plane of homogeneous modes al- 
ways constitutes an unstable focus for parameter values 
in which the domain loses stability and thus the trajec- 
tory of inhomogeneous oscillations never comes close to 
the unstable homogeneous fixed point. 




100 



FIG. 11: Thresholds for oscillatory instabilities at imposed 
current density values at which the homogeneous steady state 
becomes unstable with respect to spatial fluctuations (effec- 
tive three parameter continuation) to test condition (|l7|). 
The threshold for an oscillatory instability of the domain 
(7 hd (*o(/-0)) an d the codimension-two point (DH) analogous 
to a Turing-Hopf codimension-two bifurcation (7 DH (/i)) are 
shown as dashed and solid lines, respectively. In (a) the two 
curves are shown for the semiconductor system (double loga- 
rithmic plot) for 2d domains at the lower turning point of 
the current-voltage characteristic. The hatched re gion in- 
dicates the region in which the sufficient condition (\Vn) for 
complex spatio-temporal dynamics is fulfilled. In (bj and 
(c) 7 hd (io(/ i )) and 7 DH (/i) are shown for the electrochemi- 
cal model for the two domain bifurcations at low and high 
current densities, respectively. 



V. CONCLUSIONS 

The comparison of the two models presented in this pa- 
per allows us to identify bifurcations that exist in bistable 
systems subject to global inhibition. Apart from elec- 
trochemical and semiconductor systems such dynamics 
might be encountered in a variety of other systems, e. g. 
gas discharge devices [p8| . 

Stationary large amplitude spatial patterns called do- 
mains or filaments appear via a subcritical spatial bifur- 
cation of the uniform state and form attractors in the 
whole range of effective autocatalysis for common pa- 
rameter values in such systems. A characteristic length 
scale can be defined that facilitates quantitative compar- 
ison of the respective models. For comparable timescales 
of activator and inhibitor stable homogeneous relaxation 
oscillations can be expected. For slow dynamics of the 
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globally coupled inhibitor oscillatory instabilities of the 
domains occur, initially near the turning points of the 
current-voltage characteristic of the domain. However, 
the routes to complex spatio-temporal patterns depend 
on the local dynamics and might thus differ in each indi- 
vidual system under consideration. 

We have identified the following scenarios: If the Hopf 
bifurcation of the domain is supercritical, the system will 
display stable breathing domains. In the case of a sub- 
critical bifurcation the dynamics depends upon the fur- 
ther structure of the bifurcation diagram. If condition 
( |l7| ) is fulfilled, the onset of stable breathing or spik- 
ing modes can be expected. When inequality (^7|) is not 
fulfilled and the oscillatory instability of the domain is 
subcritical, no general statement regarding the resulting 
dynamics is possible. Either homogeneous relaxation os- 
cillations or complex spatio-temporal dynamics may re- 
sult in this case. 

We have demonstrated the above general statements 
with two models exhibiting different scenarios leading to 
stable complex spatio-temporal dynamics, thus illustrat- 
ing the general scheme. Condition (p"7|), which ensures 
that stationary or uniform modes are either unstable or 
do not exist, is fulfilled for the semiconductor system in a 
wide parameter range, but it can never be satisfied in the 
specific electrochemical model As a consequence 

the electrochemical breathing current domains always co- 
exist with homogeneous oscillations. Thus they have a 
small basin of attraction compared to the situation in the 
semiconductor model (|T^),(^4|) where no other mode is 
stable in a certain parameter range. As another conse- 
quence spiking current filaments are only present in the 
semiconductor system. Complex dynamics could only be 
found near the turning point of the current-voltage char- 
acteristic corresponding to the lower value of the imposed 
current density in both systems. 

These results emphasize the necessity to incorporate 
the spatial degree of freedom when studying electrochem- 
ical systems with negative differential resistance. Breath- 
ing current domains constitute a qualitatively new mode 
of complex spatio-temporal dynamics in electrochemical 
systems reported here for the first time. This mode may 
evolve to chaotic spatio-temporal dynamics via a period 
doubling cascade. 

It should be noted that recent experimental studies 
of the CO-oxidation on Pt-single crystal electrodes have 
shown small amplitude oscillations of the pot ential in the 
range of negative differential resistance (69|. This sys- 
tem might be an experimental illustration of the above 
results, and therefore spatially resolved measurements 
would be desirable. 



TABLE III: Typical parameter values 



k r exp(-f£V)= 



k ad = 1 
k d = 5 

N max = 1 ■ 10 
n= 1 

X= 1 
T= 300K 



10 4 
10~ 

■j^Q — 8 cm 



mol s 

3 s" 1 



C = 


20 ■ lO" 6 ^- 

V cm 


Ci = 


2-HT 6 ^-^ 

V cm^ 


Cad = 


1 ■ lO" 6 -^ 

cm 


D e = 


1 . i - 5 £si 

s 


g'= 


-1.2 • 10 5 ^-r 

mol 




1/2 



projects B4 and Bl. 



APPENDIX A: NON-DIMENSIONALIZATION 



In this section we give the transformations yielding the 
dimensionless model equations (||),(||). In physical units 
the equations read ]59| 



Co 



DL 



= i - X nFc r k r (l- < 9 >) 



• exp x 



anF 



-'DL 



-v) 



06 



kadC a d(l - 6»)exp(-aw'(0 DL ,6»)) 
-Mexp((l-aK(^ DL ^)) 



with 



1i/(0DL,0) = 



Co — Ci 
2N max k B T 



RT' 



The meaning of the numerous constants is given in |5S| 
and typical values of the constants are shown in table 
n| The model equations (g),(§) are retained via the 
transformations of the variables according to 
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and with the introduction of the parameters 



A* = 

P = 

v = 

9 = 
1 = 

10 = 



L 2 k ad c ad 
k d 



kadCad 

R 2 T(C - d) 
2N max k B an 2 F 2 
ag' 
RT 
a(LnF) 2 



cmF* 



7T 2 D«RTC ( 



c r k r exp(-x^rV) 



nFc r k r exp(- x ^V)' 



With the values given in table II] the parameters p, v 
and g are fixed to v = 0.025 @ , p = 0.5 and g = -2.4. 
They correspond to such physical values as free adsorp- 
tion sites or interaction strength. 7 depends on the well 
accessible concentration of the reacting species, which 
can be varied over several decades; [i ~ 100(L[cm]) 2 ; i 
can be set by the galvanostatic control unit and typical 
values will be of order 10 3 -10 4 . 



APPENDIX B: BIFURCATIONS AND PHASE 
PORTRAITS NEAR THE 
DH-CODIMENSION-TWO POINT 



In fig. [lj the bifurcations and phase portraits near the 
codimcnsion-two point in which the domain bifurcation 
and Hopf bifurcation of the homogeneous steady state 
meet (DH) is shown. The additional branches not shown 
in fig. ^| are a Hopf bifurcation of the unstable stationary 
domain leading to an unstable inhomogeneous limit cy- 
cle and the pitchfork bifurcation of periodic orbits that 
stabilizes the homogeneous limit cycle born in the Hopf 
bifurcation of the homogeneous steady state and which is 
the origin of another unstable inhomogeneous limit cycle 
(cf. fig. ||a). Both branches terminate in the DH. The 
respective phase portraits (insets) depict the dynamics 
schematically in a projection on the plane spanned by 
the eigenvectors of the two complex conjugate eigenval- 
ues describing the Hopf bifurcation of the homogeneous 
fixed point resp. the stationary unstable domain. The 
third direction describes the subcritical domain bifurca- 
tion (spatial mode). 
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